      implicit none
      include 'part_dist.inc'
      include 'counter.inc'
      real lbox
      character*512 input_file
      character*512 filename,position_file,dummy_file
      integer*4 imaskdim2
      integer iargc,nside
      real xc,yc,zc,side
      real*8 rcub,rat,rlen,rmtot
      integer i,j,k,ns,it0,nll0,itot,imax,ihighres
      integer ntot,ntot_old,irange,imult
      integer ix,iy,iz,iv,it,nq,nlev,nh,l,ik
      real xqdiff
      integer*8 itot_out
      integer*4 ibuffer(7)
      integer nhres,nfiles,ndiv,nrem,nlist
      character*512 pos_full_name
      integer*8 ni8,si8,ai8,ii8


      real*8 sumx,sumy,sumz,summ
      real*8 Ixx,Ixy,Ixz
      real*8     Iyy,Iyz
      real*8         Izz
      real*4 mprev
      real*4 xx0,yy0,zz0,rr0

      if (iargc().ne.1) stop 'usage: make_pl input_file_name'

      call getarg(1,input_file)
      call make_llist


      open(30,file=input_file,status='old')
      read (30,*) filename
      read (30,*) nside
      read (30,*) position_file
      read (30,*) dummy_file
      read (30,*) nfiles
      close(30)
      
      print*,'Mask file:',filename
      print*,'No. particles on a side ',nside

      open (31,file=filename,form='unformatted',status='old')
      read (31) imaskdim2,xc,yc,zc,side
      write(6,*)' imaskdim  = ',imaskdim
      write(6,*)' imaskdim2 = ',imaskdim2
      read (31) xx0,yy0,zz0,rr0
      write(6,*)' OriginalSelection_x = ',xx0
      write(6,*)' OriginalSelection_y = ',yy0
      write(6,*)' OriginalSelection_z = ',zz0
      write(6,*)' OriginalSelection_r = ',rr0
      read (31) imask
      close(31)

      if (imaskdim2.ne.imaskdim) stop 'Error: problem with imaskdim'

      print*, 'Read mask in ',xc,yc,zc,side

      do i=1,imaxv
       ival(i) = 0
      enddo

      nhres = 0

      do i=1,imaskdim
       do j=1,imaskdim
        do k=1,imaskdim
         iv = imask(k,j,i)
         ival(iv) = ival(iv)+1
        enddo
       enddo
      enddo

      ithresh(1) = 3
      ithresh(2) = 2
      ithresh(3) = 2
      ithresh(4) = 1
      ithresh(5) = 1
      ithresh(6) = 1

      if (nside.ne.imaskdim*(nside/imaskdim)) print*, 'Does not divide'

      ns = nside/imaskdim


c   Now decide how many particles in each cell of the mask

      nll(1) = ns
      nll(2) = 2
      nll(3) = 2
      nll(4) = 2
      nll(5) = 1
      nll(6) = 1
      nll(7) = 1
      do i=8,imaxv
      nll(i) = 1
      enddo


c      nll(2) = max(1,int(0.5+0.6*real(ns)))
c      nll(3) = max(1,ns/2)
c      nll(4) = max(1,ns/2)
c      it0 = ival(1)+ival(2)+ival(3)
c      it = it0
c      nll0 = max(1,nll(4)/2)
c      do i=5,imaxv
c       
c      it = it + ival(i)
c       if (it.gt.it0) then
c          nll0 = max(1,nll0/2)
c          it0 = it
c       endif
c       nll(i) = nll0
c      enddo

      print*,'nll(1) = ',nll(1)
      print*,'nll(2) = ',nll(2)
      print*,'nll(3) = ',nll(3)
      print*,'nll(4) = ',nll(4)
 

c  Now create particles

      itot = 0
      imax = imaskdim
      ihighres = 0
      ntot   = 0
      ntot_old = 0

      irange = imaskdim/ncg
      imult  = ncg 
      if (ncg*irange.ne.imaskdim) stop 'Error in values of imaskdim,ncg'

      do iz=1,irange
          print*,iz,itot
          call imaskrange
       do iy=1,irange
        do ix=1,irange


          call reset_counter

 10       continue

          if (icount(1).ne.1) goto 20

          call check_cell(ix,iy,iz,nside,ntot,itot,nhres) 

          if (ntot-ntot_old.gt.2500) then
c            print*,ntot,ntot_old,ix,iy,iz
c            stop 'Error in ntot'
          endif
          ntot_old = ntot

          goto 10

 20     continue
        enddo
       enddo
      enddo

      print*,'Number of particles used ',ntot,'/',nside**3

      print*,'Number of high res particles ',nhres

      if (ntot.ne.nside**3) stop 'Error with ntot '


c   Do a centre of mass calculation for particles in high-res cube
      sumx = 0.
      sumy = 0.
      sumz = 0.
      summ = 0.

c   Also moment of inertia calculation
      Ixx = 0.
      Ixy = 0.
      Ixz = 0.
      Iyy = 0.
      Iyz = 0.
      Izz = 0.

      do i=1,itot
        sumx = sumx + x(1,i)*ms(i)
        sumy = sumy + x(2,i)*ms(i)
        sumz = sumz + x(3,i)*ms(i)
        summ = summ + ms(i)

        Ixx  = Ixx + ms(i)*x(1,i)**2
        Ixy  = Ixy + ms(i)*x(1,i)*x(2,i)
        Ixz  = Ixz + ms(i)*x(1,i)*x(3,i)
        Iyy  = Iyy + ms(i)*x(2,i)*x(2,i)
        Iyz  = Iyz + ms(i)*x(2,i)*x(3,i)
        Izz  = Izz + ms(i)*x(3,i)*x(3,i)

      enddo

      print*,'CM: x ',sumx/summ
      print*,'CM: y ',sumy/summ
      print*,'CM: z ',sumz/summ
      print*,'CM: m ',summ

      print*,'Ixx,Ixy,Ixz: ',12.*real(Ixx/summ),12.*real(Ixy/summ),12.*real(Ixz/summ)
      print*,'Iyx,Iyy,Iyz: ',12.*real(Ixy/summ),12.*real(Iyy/summ),12.*real(Iyz/summ)
      print*,'Izx,Izy,Izz :',12.*real(Ixz/summ),12.*real(Iyz/summ),12.*real(Izz/summ)

 
c  Now add particles around the box

      lbox = 1./side
      nq = 39
      rat = real(nq+2)/real(nq)
      nlev = int(log10(lbox)/log10(rat) + 0.5)
      print*,'Number of levels: ',nlev

      lbox = rat**nlev
      print*,'Rescaling box size to: ',lbox
      nh = (nq+2)/2
      
      do l=0,nlev-1
        do k=-nh,nh
         do j=-nh,nh
          do i=-nh,nh
           ik=max(abs(i),abs(j),abs(k))
            if (ik.eq.nh) then
              rlen = rat**l
              rcub = rlen/dble(nq)
              itot = itot+1
              x(1,itot)=rcub*dble(i)
              x(2,itot)=rcub*dble(j)
              x(3,itot)=rcub*dble(k)
              ms(itot) = (rcub*dble(nside))**3
            endif
          enddo
         enddo
        enddo
       enddo


c   Do a centre of mass calculation for all particles
      sumx = 0.
      sumy = 0.
      sumz = 0.
      summ = 0.

c   Also moment of inertia calculation
      Ixx = 0.
      Ixy = 0.
      Ixz = 0.
      Iyy = 0.
      Iyz = 0.
      Izz = 0.

      do i=1,itot
        sumx = sumx + x(1,i)*ms(i)
        sumy = sumy + x(2,i)*ms(i)
        sumz = sumz + x(3,i)*ms(i)
        summ = summ + ms(i)
        Ixx  = Ixx + ms(i)*x(1,i)**2
        Ixy  = Ixy + ms(i)*x(1,i)*x(2,i)
        Ixz  = Ixz + ms(i)*x(1,i)*x(3,i)
        Iyy  = Iyy + ms(i)*x(2,i)*x(2,i)
        Iyz  = Iyz + ms(i)*x(2,i)*x(3,i)
        Izz  = Izz + ms(i)*x(3,i)*x(3,i)

      enddo

      print*,'CM: x ',sumx/summ
      print*,'CM: y ',sumy/summ
      print*,'CM: z ',sumz/summ
      print*,'CM: m ',summ

      print*,'Ixx,Ixy,Ixz: ',real(Ixx/summ),real(Ixy/summ),real(Ixz/summ)
      print*,'Iyx,Iyy,Iyz: ',real(Ixy/summ),real(Iyy/summ),real(Iyz/summ)
      print*,'Izx,Izy,Izz :',real(Ixz/summ),real(Iyz/summ),real(Izz/summ)

c  Rescale masses so that the total mass is unity
c  Rescale positions to -0.5,0.5 interval
  
      rmtot = lbox**3 * dble(nside)**3

      do i=1,itot       
      ms(i) = ms(i)/rmtot
       x(1,i) = x(1,i)/lbox
       x(2,i) = x(2,i)/lbox
       x(3,i) = x(3,i)/lbox
      enddo




c  Change the origin

c  First change xc,yc,zc to coordinate system of the parent simulation
c
c      print*,'Making HARDWIRED CHANGE to xc,yc,zc!!!!!!!!!!'
c      xc = xc + (0.51953125-0.5)
c      yc = yc + (0.28515625-0.5)
c      zc = zc + (0.23828125-0.5)


      print*,'xc,yc,zc ',xc,yc,zc


      print*,'Changing coordinate origin '
      do i=1,itot       
      ms(i) = ms(i)
       x(1,i) = x(1,i)+xc 
       x(2,i) = x(2,i)+yc
       x(3,i) = x(3,i)+zc
       do j=1,3
        if (x(j,i).ge.1.) x(j,i)=x(j,i)-1.0
        if (x(j,i).lt.0.) x(j,i)=x(j,i)+1.0
       enddo
      enddo


c      open (19,file='/data/rw2/arj/gimic/test_dir/p_part_bin',form='unformatted',
c     &    status='unknown')
c      write (19) itot
c      write (19) (x(1,i),i=1,itot)
c      write (19) (x(2,i),i=1,itot)
c      write (19) (x(3,i),i=1,itot)
c      close (19)


      print*,'Number of particles =',itot
c      print*,'Sorting by mass ... '
c      call indexx(itot,ms,ilist)

c     Establish pseudo-random ordering
      ni8 = itot
      call find_coprime(ni8,si8)
      print*,'ni8,si8 = ',ni8,si8

      do ii8 = 0,itot-1
        ai8 = mod(si8*ii8,ni8)  ! every integer 0<=ai8<ni8 generated.  
        ilist(1 + ii8) = 1 + ai8
      enddo


c      print*,'Sorting particles by x coordinate '
c
c      call indexx_x_sort(itot,x,ilist)

      print*,'Writing particle distribution to disk ...'

      itot_out = itot
      print*,itot

      ndiv = itot/nfiles
      nrem = itot - ndiv*(nfiles-1)

      do j=0,nfiles-1
      if (j.eq.nfiles-1) then
               nlist = nrem
      else
               nlist = ndiv
      endif

c      if (nfiles.eq.1) then
c        pos_full_name = position_file
c      else
       call file_name_generate(position_file,j,pos_full_name)
c      endif

    
      open (10,file=pos_full_name,form='unformatted',
     &       status='unknown')
      write (10)  nlist,itot_out,j,nfiles,ibuffer
      write (6,*) 'size(nlist)    = ',bit_size(nlist)
      write (6,*) 'size(itot_out) = ',bit_size(itot_out)
      write (6,*) 'size(j)        = ',bit_size(j)
      write (6,*) 'size(nfiles)   = ',bit_size(nfiles)
      write (6,*) 'size(ibuffer)  = ',bit_size(ibuffer)
      write (10) xx0,yy0,zz0,rr0 ! <- added for Ali's TRAPHIC
      write (6,*) xx0,yy0,zz0,rr0      
      write (6,*) nlist,itot_out,j,nfiles,ibuffer
      write (6,*) ' XXXXX ',x(1,ilist(1)),x(2,ilist(1)),x(3,ilist(1))
      write (10) (x(1,ilist(i)),i=1+j, ndiv*nfiles,nfiles),
     &           (x(1,ilist(i)),i=ndiv*nfiles+1,ndiv*(nfiles-1)+nlist)

      write (10) (x(2,ilist(i)),i=1+j, ndiv*nfiles,nfiles),
     &           (x(2,ilist(i)),i=ndiv*nfiles+1,ndiv*(nfiles-1)+nlist)

      write (10) (x(3,ilist(i)),i=1+j, ndiv*nfiles,nfiles),
     &           (x(3,ilist(i)),i=ndiv*nfiles+1,ndiv*(nfiles-1)+nlist)

      write (10) ( ms(ilist(i)),i=1+j, ndiv*nfiles,nfiles),
     &           ( ms(ilist(i)),i=ndiv*nfiles+1,ndiv*(nfiles-1)+nlist)

      close(10)

      enddo

      print*,'Particle data written to disk '

      print*,'Sorting by mass ... '
      call indexx(itot,ms,ilist)
      mprev = -200

      open (10,file=dummy_file,status='unknown')
      do i=1,itot,1001
       j = ilist(i)
       if (ms(j).ne.mprev) then
          mprev = ms(j)
          print*,i,mprev,log10(mprev)
       endif


c        if (abs(x(1,j)-0.45)<0.0025) then
c       xqdiff = abs(x(1,j)-xc)/side*20.
c       if ((xqdiff-int(xqdiff)).lt.0.1) then
          write (10,*) real(x(1,j)),real(x(2,j)),real(x(3,j)),ms(j)
c       end if
      enddo
      close(10)

      end





      subroutine highres_p(ix,iy,iz,nl,ns,nside,itot)
      include 'part_dist.inc'
      rbox = real(ns)/real(nside)
      rint = rbox/real(nl)
      rmass = (real(ns)/real(nl))**3

      do i=1,nl
       do j=1,nl
        do k=1,nl
        itot = itot+1
        x(1,itot) = real(ix-1)*rbox + (0.5+real(i-1))*rint-0.5
        x(2,itot) = real(iy-1)*rbox + (0.5+real(j-1))*rint-0.5
        x(3,itot) = real(iz-1)*rbox + (0.5+real(k-1))*rint-0.5
        ms(itot) = rmass
        enddo
       enddo
      enddo
      return
      end


      subroutine highres_glass(ihighres,nl,nside,ntot,itot,nhres)
      implicit none
      integer ihighres,nl,nside,ntot,itot,nhres
      include 'part_dist.inc'
      include 'counter.inc'
      include 'll.inc'
      real*8 xc,yc,zc
      real*8 xb,yb,zb,sz
      real*8 factor,rescale
      integer i,j,k
      


      factor =  dble(nside)/dble(ngrid)/dble(imaskdim)
      rescale = 1./factor/dble(imaskdim)
      xb = factor* (dble(ixco-1))      ! bottom-left-hand corner
      yb = factor* (dble(iyco-1))
      zb = factor* (dble(izco-1))
      if (ihighres.eq.1) then
        sz = factor*real(isize)
      else
        sz = factor*real(isize)/real(nl)
      endif

      if (ntot.eq.0) then
         print*, 'ihighres',ihighres
         print*, 'sz = ',sz
         print*,'nl = ',nl
      endif

      if (ihighres.eq.1) then
c----------------HIGH RES PARTICLES --------------------
       xc = xb + 0.5*sz
       yc = yb + 0.5*sz
       zc = zb + 0.5*sz
       call part_list(xc,yc,zc,sz)
       if (nindex.eq.0) then
           print*, 'Warning: No glass particles found!  (1)'
       endif
       do i=1,nindex
         itot = itot+1
         ntot = ntot+1
         do j=1,3
           x(j,itot) = xlist(j,i)*rescale - 0.5
         enddo
           ms(itot)  = 1.0
           nhres = nhres+1
       enddo
c--------------------------------------------------------
      else
c----------------LOW RES PARTICLES ----------------------
      do i=1,nl
       do j=1,nl
        do k=1,nl
        xc = xb + (0.5+real(i-1))*sz
        yc = yb + (0.5+real(j-1))*sz
        zc = zb + (0.5+real(k-1))*sz
        call part_list(xc,yc,zc,sz) 
        if (nindex.eq.0) then
c                print*,'factor = ',factor
c                print*,'nl = ',nl
c                print*,xc,yc,zc,sz
c                print*,i,j,k
c                print*,xb,yb,zb
c                stop 'No glass particles found! (2) ' 
c        endif 
        else            
          itot = itot+1
          x(1,itot) = xcm*rescale - 0.5
          x(2,itot) = ycm*rescale - 0.5
          x(3,itot) = zcm*rescale - 0.5

          if (ntot.eq.0) then
           print*,'xcm,ycm,zcm ',xcm,ycm,zcm
           print*,'rescale = ',rescale
           print*,itot,x(1,itot),x(2,itot),x(3,itot)
          endif
         
          ms(itot)  = nindex
          ntot = ntot + nindex
        endif 
        enddo
       enddo
      enddo
c---------------------------------------------------------
      endif

      return
      end


      subroutine vminmax(ix,iy,iz,iminval,imaxval)
      implicit none    
      include 'part_dist.inc'
      include 'counter.inc' 
      integer ix,iy,iz,iminval,imaxval
      integer i,j,k,nval
      integer ixm,iym,izm
      integer ixp,iyp,izp
      integer iix,iiy,iiz
      integer ipoint



      ixm = (ix-1)*ncg + 1
      iym = (iy-1)*ncg + 1
      izm = (iz-1)*ncg + 1
      nval = ncg

      ipoint = 2

      if (icount(2).ne.0) then


 10     continue

        call isplit(icount(ipoint),iix,iiy,iiz)
        nval = nval/2
        if (nval.eq.0) then
               print*,icount(1),icount(2),icount(3),icount(4)
               print*,ipoint
               print*,ixm,iym,izm,imask(ixm,iym,izm)
               print*,iix,iiy,iiz
               print*,ix,iy,iz
               stop 'Error in vminmax '
        endif
        ixm = ixm + nval*iix
        iym = iym + nval*iiy
        izm = izm + nval*iiz

        
        ipoint = ipoint+1
        if (icount(ipoint).ne.0) goto 10
      
      endif

      ixp = ixm+nval-1
      iyp = iym+nval-1
      izp = izm+nval-1

      iminval = 10000000
      imaxval = -1000

      do k=izm,izp
       do j=iym,iyp
        do i=ixm,ixp
            iminval = min(iminval,imask(i,j,k))
            imaxval = max(imaxval,imask(i,j,k))
        enddo
       enddo
      enddo

c  Supply cell info: (ixco,iyco,izco) give corner of cell (all -ve octant)

      ith    = ithresh(ipoint-1)
      ixco   = ixm
      iyco   = iym
      izco   = izm
      isize  = nval



      return
      end

      subroutine check_cell(ix,iy,iz,nside,ntot,itot,nhres)
      implicit none
      include 'counter.inc'
      include 'part_dist.inc'
      integer ix,iy,iz,nside,ntot,itot
      integer iminval,imaxval
      integer ihighres,nl
      integer nhres


       
      call vminmax(ix,iy,iz,iminval,imaxval)


      if (ntot.eq.0) print*,'cc: ',iminval,imaxval,ith

      if (iminval.ge.ith) then
            nl = 1
            ihighres=0
            call highres_glass(ihighres,nl,nside,ntot,itot,nhres)
            call successor
            return
      endif


      if (iminval.eq.imaxval) then          
            nl = isize*nll(iminval)
c            print*,'nl,isize,iminval ',nl,isize,iminval
            if (iminval.eq.1) then 
                 ihighres = 1
            else
                 ihighres = 0
            endif
            call highres_glass(ihighres,nl,nside,ntot,itot,nhres)
            call successor
            return
      endif



      call open_cell
      
      return
      end

      subroutine isplit(nval,iix,iiy,iiz)
      implicit none
      integer nval,iix,iiy,iiz

      if ((nval.le.0).or.(nval.gt.8)) stop 'Error in isplit'
  
      iix = (nval-1)/4
      iiy = mod(nval-1,4)/2
      iiz = mod(nval-1,2)
   
      return
      end




       subroutine reset_counter
       implicit none
       include 'counter.inc'
       integer j

          do j=1,8
              icount(j)=0
          enddo
              icount(1)=1
       return
       end

     
       subroutine successor
       implicit none
       include 'counter.inc'
       integer idigit

       idigit = 1

 10    continue

       if (icount(idigit+1).eq.0) then
 20       continue
              icount(idigit) = icount(idigit)+1
              if (icount(idigit).eq.9) then
                    icount(idigit) = 0
                    idigit = idigit-1
                    if (idigit.ne.0) then
                         goto 20
                    else
                         return
                    endif
              endif
       else
           idigit = idigit+1
           goto 10
       endif

       return 

       end


       subroutine open_cell
       include 'counter.inc'
       idigit = 1

 10    continue
  
       if (icount(idigit).ne.0) then
            idigit = idigit+1
            goto 10
       endif

       icount(idigit) = 1

       return
       end



       subroutine imaskrange
       include 'part_dist.inc'

       imin = 2147483647
       imax = 0
       do i=1,imaskdim
        do j=1,imaskdim
         do k=1,imaskdim
           imin = min(imin,imask(k,j,i))
           imax = max(imax,imask(k,j,i))
         enddo
        enddo
       enddo
       print*,'imask range:',imin,imax
       return
       end


      subroutine file_name_generate(filebase,ifile,filename)
      implicit none
      character*(*) filebase,filename
      integer ifile
      character*10 numeral
      integer lnblnk,rindex


      write (numeral,'(i10)') ifile
      filename = filebase(1:lnblnk(filebase))//'.'//
     &  numeral(1+rindex(numeral,' '):lnblnk(numeral))
      return
      end


      SUBROUTINE indexx(n,arr,indx)
      INTEGER n,indx(n),M,NSTACK
      REAL arr(n)
      PARAMETER (M=7,NSTACK=50)
      INTEGER i,indxt,ir,itemp,j,jstack,k,l,istack(NSTACK)
      REAL a
      do 11 j=1,n
        indx(j)=j
11    continue
      jstack=0
      l=1
      ir=n
1     if(ir-l.lt.M)then
        do 13 j=l+1,ir
          indxt=indx(j)
          a=arr(indxt)
          do 12 i=j-1,1,-1
            if(arr(indx(i)).le.a)goto 2
            indx(i+1)=indx(i)
12        continue
          i=0
2         indx(i+1)=indxt
13      continue
        if(jstack.eq.0)return
        ir=istack(jstack)
        l=istack(jstack-1)
        jstack=jstack-2
      else
        k=(l+ir)/2
        itemp=indx(k)
        indx(k)=indx(l+1)
        indx(l+1)=itemp
        if(arr(indx(l+1)).gt.arr(indx(ir)))then
          itemp=indx(l+1)
          indx(l+1)=indx(ir)
          indx(ir)=itemp
        endif
        if(arr(indx(l)).gt.arr(indx(ir)))then
          itemp=indx(l)
          indx(l)=indx(ir)
          indx(ir)=itemp
        endif
        if(arr(indx(l+1)).gt.arr(indx(l)))then
          itemp=indx(l+1)
          indx(l+1)=indx(l)
          indx(l)=itemp
        endif
        i=l+1
        j=ir
        indxt=indx(l)
        a=arr(indxt)
3       continue
          i=i+1
        if(arr(indx(i)).lt.a)goto 3
4       continue
          j=j-1
        if(arr(indx(j)).gt.a)goto 4
        if(j.lt.i)goto 5
        itemp=indx(i)
        indx(i)=indx(j)
        indx(j)=itemp
        goto 3
5       indx(l)=indx(j)
        indx(j)=indxt
        jstack=jstack+2
        if(jstack.gt.NSTACK)pause 'NSTACK too small in indexx'
        if(ir-i+1.ge.j-l)then
          istack(jstack)=ir
          istack(jstack-1)=i
          ir=j-1
        else
          istack(jstack)=j-1
          istack(jstack-1)=l
          l=i
        endif
      endif
      goto 1
      END


      SUBROUTINE indexx_x_sort(n,arr,indx)
      INTEGER n,indx(n),M,NSTACK
      REAL*8 arr(3,n)
      PARAMETER (M=7,NSTACK=50)
      INTEGER i,indxt,ir,itemp,j,jstack,k,l,istack(NSTACK)
      REAL*8 a
      do 11 j=1,n
        indx(j)=j
11    continue
      jstack=0
      l=1
      ir=n
1     if(ir-l.lt.M)then
        do 13 j=l+1,ir
          indxt=indx(j)
          a=arr(1,indxt)
          do 12 i=j-1,1,-1
            if(arr(1,indx(i)).le.a)goto 2
            indx(i+1)=indx(i)
12        continue
          i=0
2         indx(i+1)=indxt
13      continue
        if(jstack.eq.0)return
        ir=istack(jstack)
        l=istack(jstack-1)
        jstack=jstack-2
      else
        k=(l+ir)/2
        itemp=indx(k)
        indx(k)=indx(l+1)
        indx(l+1)=itemp
        if(arr(1,indx(l+1)).gt.arr(1,indx(ir)))then
          itemp=indx(l+1)
          indx(l+1)=indx(ir)
          indx(ir)=itemp
        endif
        if(arr(1,indx(l)).gt.arr(1,indx(ir)))then
          itemp=indx(l)
          indx(l)=indx(ir)
          indx(ir)=itemp
        endif
        if(arr(1,indx(l+1)).gt.arr(1,indx(l)))then
          itemp=indx(l+1)
          indx(l+1)=indx(l)
          indx(l)=itemp
        endif
        i=l+1
        j=ir
        indxt=indx(l)
        a=arr(1,indxt)
3       continue
          i=i+1
        if(arr(1,indx(i)).lt.a)goto 3
4       continue
          j=j-1
        if(arr(1,indx(j)).gt.a)goto 4
        if(j.lt.i)goto 5
        itemp=indx(i)
        indx(i)=indx(j)
        indx(j)=itemp
        goto 3
5       indx(l)=indx(j)
        indx(j)=indxt
        jstack=jstack+2
        if(jstack.gt.NSTACK)pause 'NSTACK too small in indexx'
        if(ir-i+1.ge.j-l)then
          istack(jstack)=ir
          istack(jstack-1)=i
          ir=j-1
        else
          istack(jstack)=j-1
          istack(jstack-1)=l
          l=i
        endif
      endif
      goto 1
      END

      subroutine part_list(xc,yc,zc,sz)
      implicit none 
      include 'll.inc'
      real*8 xc,yc,zc,sz
      real*8 xcq,ycq,zcq
      real*8 x0,y0,z0
      real*8 szqx,szqy,szqz
      integer ibox
      integer ix,iy,iz
      integer iconst
      real*8 sum


      nindex = 0
      xcm    = 0.
      ycm    = 0.
      zcm    = 0.
      sum    = 0.

      x0 = xc-0.5*sz
      y0 = yc-0.5*sz
      z0 = zc-0.5*sz

      iconst = int(sz)+2

      ibox = iconst-int(real(iconst)-sz) 


c      print*,'ibox = ',ibox

      do ix=1,ibox
       do iy=1,ibox
        do iz=1,ibox
         xcq = x0 + dble(ix-1) + 0.5*min(1., sz-dble(ix-1))
         ycq = y0 + dble(iy-1) + 0.5*min(1., sz-dble(iy-1))
         zcq = z0 + dble(iz-1) + 0.5*min(1., sz-dble(iz-1))
         szqx = min(1., sz-dble(ix-1))
         szqy = min(1., sz-dble(iy-1))
         szqz = min(1., sz-dble(iz-1))
         sum = sum + szqx*szqy*szqz
         call box_part_list(xcq,ycq,zcq,szqx,szqy,szqz)
c         print*,'nindex ',nindex,szqx*szqy*szqz
        enddo
       enddo
      enddo

      if (nindex.ne.0) then
                 xcm = xcm/real(nindex)
                 ycm = ycm/real(nindex)
                 zcm = zcm/real(nindex)
      endif

c      print*,sum
c      stop
      return
      end

      
     


      subroutine box_part_list(xc,yc,zc,szqx,szqy,szqz)
      implicit none 
      include 'll.inc'
      real*8 xc,yc,zc,szqx,szqy,szqz

      real*8 xm,ym,zm
      real*8 xp,yp,zp
      real*8 xci,yci,zci
      real*8 xdiff,ydiff,zdiff,ddif

      integer i,j,k
      integer ixcell,iycell,izcell
      integer ip,ixp,iyp,izp,iszx,iszy,iszz
      
  


      if (szqx.gt.1.) stop 'szqx>1'
      if (szqy.gt.1.) stop 'szqy>1'
      if (szqz.gt.1.) stop 'szqz>1'
      if (xc.le.0) stop 'xc<=0'
      if (yc.le.0) stop 'yc<=0'
      if (zc.le.0) stop 'zc<=0'

      xci = xc-int(xc)
      yci = yc-int(yc)
      zci = zc-int(zc)

      xm = xc-0.5*szqx
      ym = yc-0.5*szqy
      zm = zc-0.5*szqz

      if (xm.lt.0) xm=0.
      if (ym.lt.0) ym=0.
      if (zm.lt.0) zm=0.


      xp = xm - int(xm)
      yp = ym - int(ym)
      zp = zm - int(zm)

      ixp = xp*ngrid + ngrid
      iyp = yp*ngrid + ngrid
      izp = zp*ngrid + ngrid

      iszx = szqx*ngrid+1
      iszy = szqy*ngrid+1
      iszz = szqz*ngrid+1

      iszx = min(ngrid-1,iszx)  ! To avoid counting particles twice
      iszy = min(ngrid-1,iszy)
      iszz = min(ngrid-1,iszz)

c      print*,'izp,izp+iszz',izp,izp+iszz
c      print*,'iyp,iyp+iszy',iyp,iyp+iszy
c      print*,'ixp,ixp+iszx',ixp,ixp+iszx

      do k=izp,izp+iszz
       do j=iyp,iyp+iszy
         do i=ixp,ixp+iszx

            ixcell = 1+mod(i,ngrid)
            iycell = 1+mod(j,ngrid)
            izcell = 1+mod(k,ngrid)
 
            ip = igrid(ixcell,iycell,izcell)

 10         continue

            if (ip.ne.0) then
              xdiff=(xglass(1,ip)-xci)
              ydiff=(xglass(2,ip)-yci)
              zdiff=(xglass(3,ip)-zci)
              ddif=max(min(abs(xdiff),1.-abs(xdiff))/szqx,
     &              min(abs(ydiff),1.-abs(ydiff))/szqy,
     &                 min(abs(zdiff),1.-abs(zdiff))/szqz)    
              if (ddif.lt.0.5) then
                  nindex = nindex+1


                  if (xdiff.lt.-0.5*szqx) xdiff = xdiff+1.
                  if (xdiff.gt. 0.5*szqx) xdiff = xdiff-1.


                  if (ydiff.lt.-0.5*szqy) ydiff = ydiff+1.
                  if (ydiff.gt. 0.5*szqy) ydiff = ydiff-1.


                  if (zdiff.lt.-0.5*szqz) zdiff = zdiff+1.
                  if (zdiff.gt. 0.5*szqz) zdiff = zdiff-1.

                  xlist(1,nindex) = xc + xdiff 
                  xlist(2,nindex) = yc + ydiff 
                  xlist(3,nindex) = zc + zdiff 


c                  xlist(1,nindex) = xglass(1,ip)+int(xm)
c                  xlist(2,nindex) = xglass(2,ip)+int(ym)
c                  xlist(3,nindex) = xglass(3,ip)+int(zm)

c                  indlist(nindex) = ip
                  xcm = xcm + xlist(1,nindex)
                  ycm = ycm + xlist(2,nindex)
                  zcm = zcm + xlist(3,nindex)
              endif

                ip = iilist(ip)
                goto 10
            endif
         enddo
        enddo
       enddo
      

      return
      end



      subroutine make_llist
      implicit none
      include 'll.inc'
      integer i,j,k,ixp,iyp,izp

      open (10,file='glass_ascii_file',status='old')
!      open(10,file='cubic_grid',status='old')
      do i=1,nglass
       read (10,*) (xglass(j,i),j=1,3)
      enddo



      do i=1,nglass
       iilist(i) = 0
      enddo


      do i=1,ngrid
       do j=1,ngrid
        do k=1,ngrid
          igrid(k,j,i) = 0
        enddo
       enddo
      enddo

c   Make link list

      do i=1,nglass

        ixp = ngrid*xglass(1,i)+1
        iyp = ngrid*xglass(2,i)+1
        izp = ngrid*xglass(3,i)+1


        if (igrid(ixp,iyp,izp).eq.0) then
               igrid(ixp,iyp,izp) = i
        else
            j=igrid(ixp,iyp,izp)
            igrid(ixp,iyp,izp) = i
            iilist(i) = j
        endif

      enddo

      print*,'Made link list '
      end

       subroutine find_coprime(n,s)
       implicit none
       integer*8 n,s,h

       s = 0.618033989*n

 66    continue

       call hcf(n,s,h)

       if (h.ne.1) then
            s = s + 1
            goto 66
       endif

       return
       end




       subroutine hcf(s,n,h)
       implicit none
       integer*8 s,n,h
       integer*8 div
       integer*8 n1,n2

       integer ic

       n1 = s
       n2 = n

       if (s.eq.n) then
               h = 0
               return
       endif

       ic = 0

 77    continue
       ic = ic+1

       if (n1.gt.n2) then
            div = n1/n2
            n1 = n1 - div*n2
            if (n1.eq.0) then
                     h=n2
                     return
            endif
       endif

       if (n2.gt.n1) then
            div = n2/n1
            n2 = n2 - div*n1
            if (n2.eq.0) then
                       h=n1
                       return
            endif
       endif


       if (ic.gt.100) stop 'Error - too many iterations'

       goto 77

      

       end
       




